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ABSTRACT 

We report a search for fluctuations of the sky brightness toward the north 
ecliptic pole (NEP) with the Japanese infrared astronomical satellite AKARI, 
at 2.4, 3.2, and 4.1 /im. We obtained circular maps with 10' diameter field of 
view which clearly show a spatial structure on scale of a few hundred arcsec. A 
power spectrum analysis shows that there is a significant excess fluctuation at 
angular scales larger than 100" that can't be explained by zodiacal light, diffuse 
Galactic light, shot noise of faint galaxies or clustering of low redshift galaxies. 
These results are consistent with observations at 3.6 and 4.5 /im by NASA's 
Spitzer Space Telescope. The fluctuating component observed at large angular 
scales has a blue stellar spectrum which is similar to that of the spectrum of the 
excess isotropic emission observed with IRTS. A significant spatial correlation 
between wavelength bands was found, and the slopes of the linear correlations is 
consistent with the spectrum of the excess fluctuation. These findings indicate 
that the detected fluctuation could be attributed to the first stars of the universe, 
i.e. pop. Ill stars. The observed fluctuation provides an important constraints 
on the era of the first stars. 

Subject headings: cosmology: observations - diffuse radiation - early universe 
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Introduction 



Observations of the near-infrared background are believed to be very importa nt 



20021 ). A 



for investigating the formation and evolution of the first stars ( ISantos et al. 
search for the cosmic near-infrared background has been commenced using data from 
the Diffuse Infrared Background Experiment (DIRBE) /Cosmic Background Explorer 



fCOBE) and F> 


ear-in 


rared Spectrome 


ter fN] 


RS) /Infrared Telescope in Space (TRTS) 


(Hai 


Lser et al. 


1998 


; ( 


Dambresv et al. 


2001 




Arendt & Dwek 


2003; 


Matsumoto et al. 


2005 




Levenson et al. 


2007 


). The results have consistently shown that there remains 



an isotropic emission that can hardly be explained by known sources, as long as the 
same model is applied for zodiacal light. One possible origin of this excess emission is 



population HI (pop. I 



( iSalvaterra fc Ferrara 



I) stars, the first stars of the universe that caused its reionization 



2003 



Dwek et al. 



2005bl ). The observed excess emission, however, 



is much brighter than expected, and r equires a extremely high rate of star formation during 



the pop. Ill era (jMadau fc Silk 



20051 ). Furthermore, TeV-7 obser vations of distant b 



favor a low level of ne a r-infrared extraga lactic background light (IDwek et al. 



Aharonian et al. 



2006 



Raue et al. 



2005al : 



cLZclXS 



20091 ). It has been pointed out that uncertainties in 



the modelling of zodiacal light may explain the ambiguity in the re sult, since zodiacal light 



is the brightest foreground component in the near-infrared range (IDwek et al. 



2005bJ). 



In order to obtain firm results on the cosmic near-infrared background, some researchers 



looked for fluctuations ( iKashlinsky et al. 



2002 



Matsumoto et al. 



20051 ). as zodiacal light 



is very smooth. Furthermore, fluctuation measuremen ts are important for unde rstanding 



Kashlinsky et al. 



(l2005al ) reported 



the structure formation in the pop. Ill era. Recently, 
significant fluctuations at angular scales of 100" ~ 300" based on Spitzer data and 
concluded that they were caused by pop. Ill stars. However, a few studies have claimed 



that t 



2007 



le origin of t 



lis flu ct uation was the clu s tering of low redshift galaxies (jCooray et al. 



Chary et al. 



20081 ). 



Thompson et al. 



( 12007a 



) detected a significant fluctuation 
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at 1.1 /im and 1.6 fim with the Near- Infrared Camera and Multi-Object Spectrometer 
(NICMOS) /Hubble Space Telescope (HST), although their observation was restricted to 
angular scales smaller than 100". However, the y claimed that the observ ed fluctuation 



2007bJ). In this paper, 



could be due to faint galaxies at redshift z < 8 ([Thompson et al. 
we present new observations of fluctuations in sky brightness using the Japanese infrared 
astronomical satellite AKARI. AKARI has a cold shutter to measure dark current precisely 
as well as a short wavelength band at 2.4 /im, which provided new observational evidence 
on the study of the sky fluctuation. 

This paper is organized as follows. In §2, we provide detailed information on the 
observation and data reduction procedure. In §3, we describe the fluctuation analysis of the 
near-infrared sky, and we discuss the possible origins of the fluctuation in §4. Discussions 
on the implication of our data in view of the first stars are given in §5. The final section 
summaries our results. 



Observations and data reduction 



ed 68 cm apertur e 



AKARI is the first Japanese infrared astronomy satellite with a coo 
telesc ope, launched into a Sun- synchronous orbit on February 22, 2006 dMurakami et al. 
2007). The Infrared Camera (IRC) is one of the instruments onboard AKARI (jOnaka et al 
20071 ). and has a frame of 10' x 10' with a pixel scale of 1.5". The IRC performed a 



large area survey around the north ecliptic pole (NEP) (ILee et al. 



a fixe d position in the sky near the NEP to evaluate its performance (IWada et al. 



2009). and monito red 



20071 ). The observed region of the sky was named NEP monitor field (02000 =17h55m24s, 
^2000=66:37:32). Pointed observations were carried out twice a month, and three sets of 
imaging observations were performed for all wavelength bands of the IRC during one 
pointed observation. Dithering of ~ 10" was applied between these three sets of imaging 
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observations. Since maneuvering of the satellite sometimes commenced before completing 
the third set of imaging observations, the last image at 4.1 fim was usually useless. The 
integration time for one image was 44.41 seconds, and the data were retrieved with Fowler 
4 samplings. Systematic errors due to uncertainty in the the absolute calibration were 2.8, 
2.5, and 3.4 %, which were small compared to other e rrors. The details of o bservation mode 



and calibration are shown in the IRC User's Manual (ILorente et al. 



2008|). 



In order to detect the fluctuations in sky brightness, we used near-infrared images 
taken in 14 pointed observations between September 2006 and March 2007, during which 
no radiation from the earth entered the baffle and stray light was negligible. Incidence of 
earthshine into the telescope tube was detected from the temperature sensor on the baffle 
and confirmed by the IRC simultaneously as a spurious stray radiation. Furthermore, we 
did not use images when the attitude control was less stable than usual, as well as those 
images suffered from a large number of cosmic ray hits. After all, we chose 40, 39, and 
24 images from 2.4, 3.2, and 4.1 /j,m bands, respectively, for the analysis, (hereafter we 
list data at these three bands in order of wavelength). Total integration times were 29.6, 
28.9, and 17.8 minutes fo r each band, which are rather shorter than Spitzer observations 



teashlinsky et al. 



2005al ). IRC/AKARI, however, has a larger field of view, 10' than 
IRAC/Spitzer, 5.2'. Since pixel throughputs are almost same, IRAC/Spitzer needs ~ 4 
times more observing time than IRC/AKARI to observe the same area of the sky with the 
same S/N ratio. 



We analyzed the sky image following the IRC User's Manual (ILorente et al 



2008|) 



with some improvements. We estimated the dark current of the image field from that of the 



maske d region on the array by applying the new method developed by 



Tsumura fc Wada 



( 120 111 ). This method provides a more reliable estimate of the dark current for those frames 



suffering from an after effect due to strong cosmic ray hits occurring in the South Atlantic 
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Anomaly region. We also constructed better flat fields using all images obtained during the 
period as used for the observation of the NEP monitor field. In total, 774, 771, and 516 
images were stacked, and significantly better flat fielding was achieved compared to the 
ordinary pipeline. For the correction of the aspect ratio, we used the "nearest" option as 
interpola tion type in the MAGNIFY task of the Image Reduction and Analysis Facility 



(IRAF) ( jTody Ill986l ) to avoid a stripe-like artifact. 



Finally we stacked the images and obtained the maps of the NEP monitor field. Due 
to constraints on attitude control, the position angles of the images were rotated by 1° 
per a day, and by 200° during the observation period, resulting in circular images with 10' 
diameters. This improves the quality of the stacked images compared to usual stacking, 
since flat field errors, fluctuations in sensitivity between pixels, residuals in dark frame 
subtraction, fluctuation of the zodiacal light (if any), etc. are reduced, except for the central 
part of the image and any axially symmetric structure. We also obtained dark images with 
the cold shutter closed. A dark image corresponding to each sky image was retrieved before 
or after the pointed observation. As a result, we obtained equal numbers of sky and dark 
images for each band. For the dark images, we applied the same procedure as for the sky 
images and stacked them to obtain dark maps for the three wavelength bands. 

In order to obtain the fluctuations in sky brightness, we masked pixels whose signal 
exceeded 2 a level. We repeated this 10 times after which additional clipping did not show 
significant changes. However, there remained bright pixels around the masked region due 
to the extended point spread function (PSF). In order to subtract this effect, we searched 
for foreground sources with a detection threshold of 2 a in the unmasked stacked images, 
and removed their contributions using the DAOPHOT package in the IRAF. At first we 
carefully modeled the PSF based on the beam profiles of the bright point sources in the 
unmasked stacked images. Using the clearly identified point sources (sources well fitted by 
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the PSF), we estimated the contribution of sources in the outer region using the model PSF 
and subtracted them from the sky images. We were left with some point sources (sources 
rejected by the PSF fit) that can not be fitted with the PSF due to faintness, overlapping of 
sources or other reasons. For these sources, we simply masked 8 neighboring pixels around 
their centers . We identified extended objects using ground-based data for the NEP region 



( jlmai et al. 



20071 ) which have a better spatial resolution. We created convolved images 
using AKARI's PSF and subtracted them from the sky images. For these PSF subtracted 
images, we masked the same region in the 2 a clipping image. In order to minimize 
contributions from PSFs, we further masked a layer of one pixel around the masked regions 
in which the pixel of 2 o brighter level exists. For the dark maps, we masked the same 
region as with the sky maps and repeated 2 a clipping 10 times. Table 1 shows some 
parameters used and obtained in this analysis. Percentages of remaining pixels are ~ 47% 
for all wavelength bands, and average sky brightness are 114, 73, and 105 nWm _2 sr _1 . As 
for the sky brightness, systematic errors are not taken into account in order to show the 
statistical significance clearly. 

Figure 1 shows the final fluctuation maps. The upper and lower panels indicate the 
sky and dark maps, respectively. The sky maps clearly show structure, while the dark 
maps look like random noise. As it is well known, a major emission component in the 
near-infrared sky is zodiacal light. We have to keep in mind the fact that the contribution 
of galactic stars is negligible, but there remains integrated light of galaxies fainter than the 
limiting magnitudes, and isotropic emission. 



3. Fluctuation spectrum 



Fluctuation analys is was first performed by iKashlinsky et al. I (120021) an d 

Odenwald et al. I ( 120031 ) using 2MASS data, and followed by IKashlinsky et al. I (l2005al . 
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2007ah : lArendt et al. I (12010( 1 for Spitzer data and by iThompson et al. I fl2007al Jbl) for HST 
data. In order to examine the angular dependence of the fluctuation observed with AKARI, 
we performed a power spectrum analysis for both the sky and the dark maps (Figure 1), 
following the works of previous authors. A two-dimensional Fourier transformation was 
computed for the fluctuation field 8F(x) as a function of the two dimensional angular 
wavenumber, q, and the two dimensional coordinate vector of the image, x. 



F(q) = J 6F(x) exp ( - ix - q) d 2 x 



The two-dimensional power spectrum and fluctuation spectrum are given by 
P 2 (q) = |-^(q)| 2 , and [g 2 P2(q)/(27r)] 1 / 2 . Figure 2 shows the amplitude maps of the 
two-dimensional fluctuation spectra in Fourier space. The maps are symmetric but faint 
annular patterns appear in 2.4 and 3.2 /im maps. These patterns correspond to small peaks 
in the fluctuation spectra (Figure 3) at ~ 5" and are probably artifacts caused by the 
masking pattern, since the same features appear in both sky and dark images. However, 
they do not affect the conclusion of this paper. The one dimensional fluctuation spectrum, 
W'P'il (27I")] 1 / 2 , is obtained by taking ensemble averages in the upper halves of the concentric 
rings of Figure 2 with equal spacing in the angular wavenumber, q. 

Figure 3 shows the fluctuation spectra. The upper panel shows the fluctuation spectra 
for both the sky (filled circles) and dark (open triangles) maps. Statistical errors, given 
by standard deviations divided by the square root of the number of the data points, are 
plotted. The fluctuation spectra of the sky maps are significantly larger than those of the 
dark maps. 

We subtracted the fluctuation spectra of the dark maps from that of the sky maps 
in quadrature, and show the result as filled circles in the lower panel of Figure 3. The 
solid lines indicate the fluctuation spectra due to the shot noise of faint galaxies that were 
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estimated by simulation. We first obtained galaxy counts for the stacked images before 
masking and extended them to the faint end assuming a slope of 0.23 for logarithmic counts 



( IMaihara et al. 



20011 ). We constructed fluctuation maps by randomly distributing the 
galaxies fainter than the limiting magnitudes, and performed the power spectrum analysis. 
The fluctuation spectra of the shot noise of faint galaxies fit well the observed power spectra 
with limiting magnitudes (AB) of 22.9, 23.2, and 23.8. The decline at small angle is caused 
by the PSF, while the absolute value depends on limiting magnitude and galaxy counts. 
Figure 3 shows that a significant residual fluctuation above the shot noise level remains at 
large angular scales. Table 2 shows the numerical values of fluctuation of the observed sky, 
fluctuation due to shot noise by faint galaxies and excess fluctuation at angles larger than 
100". 

We performed a subset analysis to confirm the celestial origin of the observed structure. 
We divided our set of images into two subsets corresponding to two cases. In case A, each 
subset is taken as a sequence alternating in time, while in case B, the subsets correspond to 
the earlier and later halves of our set. For both cases, we obtained two stacked images, Fl 
and F2, by applying the same procedure as the one described previously in this section. The 
fluctuation spectra of the difference between these two stacked images are shown in Figure 
4. Case A and B are shown as squares and asterisks, respectively, while the fluctuation 
spectra of the dark maps are shown by triangles. The results for both subsets are consistent 
with those of the stacked dark maps. This indicates that the observed structure is indeed 
present in the original images and is of celestial origin. 

We also examined the impact of masking on the fluctuation spectra, since ~ 53 % of all 
pixels are masked. We constructed a common mask that includes all pixels masked in any 
of three wavelength bands. The fraction of remaining pixels in images with the common 
mask applied is ~ 32 %. We again performed the analysis for these images, and compared 
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the obtained fluctuation spectra with the original ones shown in Figure 5. The fluctuation 
spectra of the images using the common mask have slightly lower than the original ones, 
but are consistent with the original ones within the error bars. This result assures that the 
masks adopted do not cause significant artifacts and that the detected fluctuation is a real 
structure in the sky. 

The upper panel of Figure 6 shows smoothed sky maps that were obtained by taking 
averages within circles of 50" diameter centered at each pixel. The smoothed maps clearly 
show the existence of a significant extended spatial structure with an angular scale of a 
few hundred arcsec for all wavelength bands, and the overall pattern is very similar. This 
is consistent with the excess fluctuation at large angular scales shown in Figure 3. The 
smoothed dark maps in the lower panel of Figure 6 do not show any structure besides 
random noise, thus they reconfirm the celestial origin of the observed large scale structure. 

In order to examine the fluctuation quantitatively, we took the averages of the excess 
fluctuation at angular scales between 100" and 350". Figure 7 shows the spectrum of 
the fluctuating component (filled circles) with 1 a errors. The characteristic feature 
of the fluctuating component is a blue-star-like spectrum, which is consistent with a 
Rayleigh- Jeans spectrum, ~ A~ 3 (straight line in Figure 7). We note that the observed 



spectrum of the fluctuating component is simi 



emission observed by IR TS (IMatsumoto et al 



at 3.6 /im and 4.5 /im (IKashlinsky et al- 



ar to the spectrum of the isotropic excess 



2005a . 



2005). Open squares show the Spitzer results 



2007al ). which are basically consistent 



with ours. We did not plot the Spitzer results at longer wavelength bands, since the data 
were scattered an d appeared less rel i able. A t the wavelength bands at 1.1 /im and 1.6 /im 



(NICMOS/HST), 



Thompson et al. 



(12007br ) detected fluctuation of 0.3 and 0.4 nWm^sr" 1 , 
respectively, at an angular scale of 85". These are lower than the values extrapolated from 
the AKARI data, assuming the Rayleigh- Jeans spectrum, ~ A~ 3 . However, the angular 
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scale for NICMOS/HST is too small to have detected the large-scale structure shown in the 
smoothed sky maps (Figure 6). Therefore, we consider it is inappropriate to compare the 
NICMOS/HST results directly with those of AKARI and Spitzer. 

We performed a two-points correlation (autocorrelation) analysis, Cf(Q) = 
(5F(x + 9)5F(-x)), both for the sky and dark maps. In order to save computation time, 
we reduced the angular resolution by adopting the average values groups of 2 x 2 pixels. 
The results, C 2 a(0), C 3 .2(#), and C 4 .i(#) are shown in Figure 8 as a function of 9 in units of 
arcsec. The data were in groups of 4 pixels, corresponding to an angular resolution of 6". 
Error bars indicate statistical errors, i.e. standard deviation divided by the square root of 
the number of pairs. It must be noted that the sampling at angles larger than 300" is not 
complete, that is, central parts of images are not counted. Figure 8 clearly shows excess 
low frequency components which are common for all three wavelength bands. The low 
frequency components correspond to the excess fluctuations at large angles, as shown in 
Figure 3. 

We made a point-to-point correlation study between the wavelength bands using 
independent data in smoothed maps. The result is shown in Figure 9. The correlation 
between 3.2 and 2.4 /im is fairly strong (correlation coefficient ~ 0.73), while that between 
4.1 and 2.4 fim is somewhat weak (correlation coefficient ~ 0.42). This is probably because 
shot noise due to faint galaxies becomes important at 4.1 /im. Colors of the correlating 
component (slope of Figure 9) normalized to be 1.0 at 2.4 /im (right ordinate) are also 
plotted in Figure 7 as open circles, which is similar to that of large scale fluctuation (filled 
circles) . 

We also obtained the cross correlation function, C_f®g(#) = (SF(x + 9)5G(x)) between 
the maps taken at different wavelengths. C2Am.2{9) and 62.4^4.1 (#) are shown in Figure 10. 
Degradation of angular resolution and binning of the data are the same as for the two-point 
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correlation analysis. Error bars indicate statistical errors, i.e. standard deviation divided 
by the square root of the number of pairs. The cross correlation functions in Figure 10 
are very similar to the 2.4 /zm two-points correlation, 6*2.4(6'). This is consistent with a 
point-to-point correlation analysis (Figure 9), and the ratios of cross correlation functions 
to 62.4(0) are similar to the slopes in Figure 9. The result of our cross correlation analysis 
implies that the fluctuation patterns are basically the same for all three wavelength bands. 



4. Origin of the excess fluctuation 



Being the main component of sky brightness, zodiacal light is a candidate for the origin 
of the fluctuation. The subset analysis provides evidence that zodiacal light is not the 
source of the observed fluctuation. Particularly in case B, the interplanetary dust observed 
during the first half period is almost independent from the one seen in the latter half period, 
since the interplanetary dust in the column along the line of sight has different orbits and 
eccentricities. If the observed fluctuation were caused by zodiacal light, the fluctuation 
spectra of the differences between these two subsets should be higher than those of the dark 
images. 



The detected fluctuation 



IRAS flVrtilek fc Hauser 



evels are also higher than tho se expected for zodiacal ligh t 



1995h . COBE (IKelsall et al. 



19981 ) and ISO (lAbraham et al. 



19971 ) attempted to detect the fluctuation of zodiacal emission on scales from arcmin to 



degrees. The most stri ngent upper limit of 0.2 % of the sky brightness was obtained by 



Abraham et al. 



(119971 ) who measured the fluctuation o : 



45' x 45' field with 3' pixel scale. Recently, 



Pyo et al. 



zodi acal emission at 25 /im for a 



(120 111 ) analyzed the smoothness of 



the mid-infrared sky using the same data set as the one we used in this paper. They applied 
the same power spectrum analysis for the mid infrared bands, but did not detect any excess 
fluctuation over the random shot noise caused mainly by photon noise. We, therefore, 
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assume the lowest fluctuation level at the 18 fj,m band to be an upper limit of fluctuations of 
the sky at an angle larger than 150", that is, ~ 0.02% of the sky brightness. Although this 
upper limit includes fluctuation due to zodiacal emission, diffuse Galactic light, integrated 
light of galaxies and cosmic background, we simply assume that the fluctuation of the 
zodiacal emission is lower than ~ 0.02% of the absolute sky brightness. 

Fluctuation of zodiacal light is similar to that of z odiacal emission. In terplanetary 



dust clouds are thought to consist of three components (IKelsall et al. 



19981 ) . The first is a 



smooth cloud, which is broadly distribut ed in the solar system . The second is dust bands 



composed of asteroidal collisional debris (jNesvorny et al. 



20031 ). The third is a circumsolar 



ring which is composed of dust resonantly trapped into orbit near 1 AU. We calculated 
the contributi o ns of these components to the integrated brightness using the model by 



Kelsall et al. 



(1998). Figure 11 shows the volume emissivity of zodiacal light at 2.2 jum 
(scattering part, solid line) and that of zodiacal emission at 12 /j,m (thermal part, dashed 
line) along the line of sight towards the NEP. The volume emissivity is normalized such 
that the total emissivity at the origin is 1.0. Contributions of the dust bands (bottom), the 
circumsolar ring (middle) and the total emissivity (top) are shown separately. We used the 
model at the epoch 2006-09-15 when the contribution by the circumsolar ring reaches its 
maxumum. Figure 11 shows that the smooth cloud is a main contributor for both zodiacal 
light and emission, while that of the dust bands is very low although they are conspicuous 
at low ecliptic latitudes. Figure 11 clearly shows that the spatial distribution of the volume 
emissivity is almost the same for scattered light (zodiacal light) and thermal emission 
(zodiacal emission). This is because scattering phase function and temperature of the dust 
in the column do not change much in this direction. Although there is a small deviation 
for the total emissivity at the distance of 1 AU and beyond, and a systematic difference 
between scattering and thermal part for the circumsolar ring, their contribution to the total 
surface brightness is only a few percent. We conclude that the difference of the fluctuation 
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between zodiacal light and zodiacal emission is at most a few percent, even if there exists a 
fluctuation of the dust density along the line of sight. 

More quantitatively, we estimated fluctuation levels to be expected due to zodiacal 
light by constraining the upper limit of the fluctuation of zodiacal light to 0.02 %. The sky 
brightnesses corresponding to these upper limits are 0.023, 0.015, and 0.021 nWm _2 sr _1 . 
The expected fluctuation of zodiacal light could be reduced further during the stacking 
procedure, since AKARI moved along the earth orbit and IRC observed different regions of 
the solar system during the observation period. The observed fluctuations at large angular 
scales, 0.25, 0.11, and 0.033 nWm~ 2 sr -1 (see Figure 7), are significantly larger than those 
expected for zodiacal light. 

The spectrum of fluctuation observed with AKARI shows a very blue color up to 
4.1 fim (Figure 7), while IRTS observations d elineated that zodiaca l emission becomes 



prominent at wavelengths longer than 3.5 /iin (jOotsubo et al 



19981 ). This difference also 



indicates that the fluctuation observed by AKARI is not due to zodiacal light. 

Based on these arguments, we conclude that the contribution of zodiacal light to the 
fluctuation is negligible. 

The second possible origin of the fluctuation is diffuse Galactic light (DGL), i.e. 
scattered starlight and thermal emission of the interstellar dust. Since DGL is closely 
related to the dust column density, we examined its correlation with the far-infrared map. 
The Far-Infrared Surveyor (FIS) / AKARI detected a clear cirrus feature toward NEP at 90 



/im with ~ 30" angular resolution (IMatsuura et al. 



2010). Figure 12 shows the correlation 



between far-infrared and near-infrared maps with the angular resolution of the AKARI 
data degraded to 30". No significant correlation was found, indicating that DGL is not the 
source of the observed near-infrared fluctuation. 
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The third candidate is the clustering of faint galaxies. A few authors (jCooray et al. 



2007 



Chary et al. 



20081 ) have claimed that the fluctuation observed with Spitzer is due 



to the clustering of low redshift galaxies. Th e blue color observed by AKARI is not 



consistent with the proposal by 



Chary et al. 



component originates in red dwarf galaxies. 



(200 81) which st a tes tha t the fluctuating 



Kashlinsky et al. 



( I2007bl ) carefully examined 
the correlation between the map of these galaxies and those observed with Spitzer and 



conclu ded that red dwarf galaxies are not responsible for the Spitzer result. 



Sullivan et al. 



(120071 ) examined the contribution of the clustering of faint galaxies to the sky fluctuation 
using the catalog of galaxies. Based on their work, we estimated that the fluctuation at 2.4 
/jin caused by clustering of galaxies fainter than the limiting magnitude of this study at 
multipole, I ~ 1,000 (~ 600"), to be ~ 0.03 nWm~ 2 sr~ 1 . This is s ignificantly lowe r than 



the fluctuation detected at 2.4 //m. As for the L band, Figure 8 in 



Sullivan et al. 



(120071 ) 



shows that the fluctuation due to clustering is fairly weak compared with the one observed 



by Spitzer (IKashlinsky et al. 



2005a 



2007al ). These results indicate that the clustering of 



low redshift faint galaxies is not likely to be the source of the observed fluctuation. 



Contribution of faint galaxies to the 


Dackground is studied 


by several aut 


lors 


Totani & Yoshii 


2000; 


Totani et al. 


2001 




Maihara et al. 


2001; 


Keenan et al. 


2010) 



According to the recent observation by 



Keenan et al. 



(120101 ) , the integrated light of galaxies 



in K band amounts to 10.0 ± 0.8 nWm 2 sr . Galaxies of ~ 18 mag are mostly responsible 
for the integrated brightness, and the galaxies fainter than that contribute even less. The 



integrated brightness at 2.4 /xm due to galaxies 



mag is estimated to be at most 0.7 nWm 2 sr 1 (IKeenan et al 



ainter than the limi ting magnitude of 22.9 



2010j). Ob served fluctuatio n 



at lar ge angular scales, 0.25 nWm 2 sr 1 , is a fraction of this brightness. 



Sullivan et al. 



(120071 ) estimated the fluctuation due to clustering of galaxies with K magnitudes between 



17 and 19 mag. The fluctuation o btained was 



of these galaxies is 4 nWm 2 sr 1 (IKeenan et al 



nWm 2 sr x , while the integrated light 



2010 ), that is, the fluctuation is 1/40 
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of the integrated light of galaxies. Furthermore, the fainter galaxies contribute less to 
the clustering. Compared with this result, the fluctuation observed with AKARI is fairly 
l arge. The situation is m ore conspicuous for the case of NICMOS/HST observations 



( iThompson et al. 



2007bl ). since their limiting magnitudes, 28.5 mag, are much fainter than 



the one of the AKARI observation. T 



re integrated light in H band due to galaxies fainter 



than 28.5 mag is ~ 0.04 nWm 2 sr 1 (IKeenan et al. 



20101 ) . while the observed fluctuation 



is 0.4 nWm _2 sr _1 . The observed fluctuation is ten times brighter than the integrated 
light of galaxies, which indicates that the origin of the observed fluctuation is not the low 
redshift galaxies. These considerations imply that a simple extrapolation of galaxy counts 
hardly explains fluctuation observed, and that there must exist a new component of infrared 
sources with magnitudes fainter than the current limit. 

Accreting black holes could be another candidate if they are abundant in the early 
universe, since they show blue spectra. However, the bulk of the light will be in much 
shorter wavelengths like X-rays, since the accretion disks must be very hot. 

To summarize, we found that the observed near-infrared fluctuation at large angular 
scales can hardly be explained by known sources. We now regard that the most probable 
origin of the fluctuation is pop. Ill stars, the first stars of the universe. 



5. Discussions 

AKARI observations have shown that there existed a large scale structure even at the 
pop. Ill era. Since the fluctuation power is nearly constant at scales of a few hundred 
arcsec, the large scale structure at the pop. Ill era could be more extended. The co-moving 
distance for 10' at z ~10 is ~ 30 Mpc which is a scale similar to the local super cluster in 
the present universe. Theoretical studies of sky fluctuation due to the first stars have been 
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2010h . 



predi cted by several authors (jCooray et al. 112004 ; iKashlinsky et al. 112004 ; iFernandez et al. 



Cooray et al. 



( 120041 ) presented the first theoretical study on fluctuations caused by 
pop. Ill stars. They assumed biased star formation which traces the density distribution 
of the dark matter. Although the absolute level of the fluctuation spectrum has a wide 
range depending on the physical parameters, the overall shapes have almost the same 
convex feature with a flat peak at / ~1,000 (~ 10') and a turn-over towards large angles. 
We regard this flat peak as the most significant excess fluctuation observed by AKARI 
which is mainly due to the pop. Ill stars. The fluctuating po wer observed with AK ARI is 



Cooray et al. 



200J). The 



marginally consistent with a pessimistic estimate (Figure 3 in 
wavelength dependence of fluctuation power of the model shows fairly blue color which is 
basically consistent with AKARI observation. The fluctuation observed with AKARI could 
be interpreted within the framework of this model. However, detection of the turn-over at 
angular scale larger than 10' is essential. 



Recently, IFernandez et al. I (120101 ) made a numerical simulation of the cosmic 
near-infrared background caused by early populations (pop. II and III stars), assuming 
a halo mass of 2 x 10 9 M & . They calculated both amplitude and fluctuation, taking 
into account stellar emission and emission from ionized gas surrounding stars. As for the 
fluctuation spectrum, they predicted a monotonic decrease of the fluctuation towards large 
angles, which is not consistent with our result. This discrepancy requires some kind of 
modification of the theory. However, their result provides valuable information on the 
physical conditions of early populations which are independent from halo mass. 



Fernandez et al. 



(2010) 



A spectrum of the sky brightness at 2 ~ 4 /im obtained by 
shows a blue color, and the main emission component is attributed to stars. Furthermore 
the ratio of the fluctuating power to the absolute sky brightness, 51/1, does not depend 
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on the wavelength. This justifies to regard the spectral shape of absolute 
same as that of the fluctuation. In Figure 7, the spectrum of Figure 20 of 



brightness to be 



Fernandez et al. 



( 120101 ) is plotted as a dotted line. Their spectrum is a little redder than, but is qualitatively 
consistent with, the AKARI observation. A good correlation for the fluctuations between 
the wavelength bands observed with AKARI supports the idea that the emission source 
is stellar emission from early populations, and that t he nearest stars provide the largest 



Fernandez et al. 



(120 lOf ) is ~ 0.02 



contribution to the sky brightness. Sl/l estimated by 
at a few hundred arcsec for the largest fluctuation. The observed fluctuation at 2.4 
/im, ~ 0.25 nWm _2 sr~ 1 , renders the absolute sky brightness to be brighter than ~ 13 
nWm~ 2 sr~ 1 , which is con sistent with the excess brightness observed with IRTS and COBE 



( jMatsumoto et al. 



20051 ) . 



Al 



2005b 



Cooray et al. 



2004 : 



(Santos et al. 


200^ 


J; Sa 


vaterra & Ferrara 


2003; 


Dwek et al. 


Fernandez et al. 


2010 


) predict a clear peak due to the 



redshifted Lya from the ionized gas surrounding first stars. The wavelength of this peak 
depends on the redshift of the nearest first stars. Since there is no signature of redshifted 
Lya in AKARI data, it must appear below 2 /im, that is, the redshift of the nearest first 
stars is less than 15. Detection of the redshifted Lya both in the amplitude and the 
fluctuation will be a key issue to delineate the epoch of formation of the first stars. 

The fluctuation observed with AKARI provides an evidence for light from the first 
stars. However, more observational results are required in order to investigate the physical 
conditions and the star formation history of the first stars. In a forthcoming paper, we 
will present a fluctuation analysis at larger angles than this paper based on the AKARI 
NEP wide field survey. The sounding rocket experiment, Cosmic Infrared Background 
Experiment (CIBER), will provide the spectrum of the sky brightness aroun d 1 /im and 



the fluctuation at 1° scale at the 0.8 /xm and 1.6 /im band (IBock et al 



20061 ). The Korean 
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infrared satellite, Multi Purpose Infrared Imaging System (MIRIS) (lHan et al. 1120101 ) will 
observe the fluctuations at even larger angular scale (~ 10°). These space observations will 
certainly improve our understanding of the first stars of the universe. 



6. Summary 

We performed a fluctuation analysis for a sky region with a diameter of 10' towards 
the NEP at 2.4, 3.2, and 4.1 jum based on observations with AKARI. A power spectrum 
analysis indicates that there exists a significant excess fluctuation at angular scale larger 
than 100" that can't be explained by zodiacal light, diffuse Galactic light, shot noise from 
faint galaxies or clustering of low redshift galaxies. The fluctuating component at large 
angles has a blue stellar spectrum, and shows a significant spatial correlation between 
adjacent wavelength bands. 



2007bl ) 



The fluctuation observed with AKARI and NICMOS/HST ([Thompson et al. 
is considerably higher than the one expected from galaxy counts. This implies that a simple 
extrapolation of galaxy counts hardly explains the observed fluctuation, and that there 
must be a new component of infrared sources at the faint end. 

The most probable origin of excess fluctuation is pop. Ill stars, the first stars of the 
universe. The detec ted fluctuation can be explained with biased star formation which 



traces dark matter (jCooray et al. 



20041 ) . A blue color of the observed fluctuation and a 



significant spatial correlation between wa velength bands imply th at the source of emission 



is stellar emission from early populations ([Fernandez et al. 



20101). 



The angular scale of 10' at z ~10 corresponds to a comoving distance of 30 Mpc, a 
scale similar to local super cluster in the present universe. 



The detection of a turn-over of fluctuation spectra at large angles and the spectral peak 
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of redshifted Lja at wavelengths shorter than 2 /iin is a key issue for future observations in 
order to improve the investigation of the pop. Ill star formation. 

ACKNOWLEDGMENTS 

We thank the IRC/AKARI team for their encouragement and helpful discussions. 
Thanks also to Professors A. Cooray, E. Komatsu and A. Ferrara for valuable discussions 
and comments. This work was supported by JSPS Grant-in-aid (18204018). HML was 
supported by NRF grant No. 2006-341-C00018. 



- 21 - 

REFERENCES 

Abraham, P., Leinert, C, & Lemke, D. 1997, A&A, 328, 702 
Aharonian, F., et al. 2006, Nature, 440, 1018 
Arendt, R.G., & Dwek, E. 2003, ApJ, 585, 305 

Arendt, R.G., Kashlinsky, A., Moseley, S. H., k Mather, J. 2010, ApJS, 186, 10 
Bock, J., et al. 2006, New A Rev., 50, 215 

Cambresy, L., Reach, W. T., Beichman, C. A., & Jarrett, T. H. 2001, ApJ, 555, 563 
Chary, R., Cooray, A., & Sullivan, I. 2008, ApJ, 681, 53 
Cooray, A., et al. 2004, ApJ, 606, 611 
Cooray, A., et al. 2007, ApJ, 659, L91 

Dwek, E., Krennrich, F. & Arendt, R. G. 2005, ApJ, 634, 155 

Dwek, E., Arendt, R. G., & Krennrich, F. 2005, ApJ, 635, 784 

Fernandez, E., Komatsu, E., Iliev, I. T., & Shapiro, P. R. 2010, ApJ, 710, 1089 

Han, W., et al. 2010, Proc. SPIE, 7731, 55 

Hauser, M. G., et al. 1998, ApJ, 508, 25 

Imai, K., et al. 2007, AJ, 133, 2418 

Kashlinsky, A., et al. 2002, ApJ, 579, L53 

Kashlinsky, A., et al. 2004, ApJ, 608, 1 

Kashlinsky, A., Arendt, R. G., Mather, J., & Moseley, S. H. 2005a, Nature, 438, 45 



-22 - 

Kashlinsky, A. 2005b, Phys. Rep., 409, 361 

Kashlinsky, A., Arendt, R. G., Mather, J., & Moseley, S. H. 2007a, ApJ, 654, L5 
Kashlinsky, A., Arendt, R. G., Mather, J., & Moseley, S. H. 2007b, ApJ, 666, LI 
Keenan, R. C, Barger A. J., Cowie L. L., & Wang, W. -H. 2010, ApJ, 723, 40 
Kelsall, T., et al. 1998, ApJ, 508, 44 
Lee, H. M., et al. 2009, PASJ, 61, 375 

Levenson, L. R., Wright, E. L., & Johnson, B. D. 2007, ApJ, 666, 34 

Lorente, R., et al. 2008, AKARI IRC Data User Manual (version 1.4), 

http:/ /www. ir.isas.jaxa.jp/ASTRO-F/Observation/DataReduction/IRC/ 

Madau, P., & Silk, J. 2005, MNRAS, 359, L37 

Maihara, T., et al. 2001, PASJ, 53, 25 

Matsumoto, T., et al. 2005, ApJ, 626, 31 

Matsuura, S., et al. submitted to ApJf arXiv: 1002 .3674) 

Murakami, H., et al. 2007, PASJ, 59, S369 

Nesvorny, D., Bottke, W. F., Levison, H. F., & Dones, L. 2003, ApJ, 591, 486 

Odenwald, S., Kashlinsky, A., Mather, J., Skrutskie, M. F., & Cutri, R. M. 2003, ApJ, 583, 
535 

Onaka, T., et al. 2007, PASJ, 59, S401 

Ootsubo, T., et al. 1998, Earth, Planets and Space, 50, 507 



-23 - 

Pyo, J., et al. 2011, submitted to ApJ 

Raue, M., Kneiske, T., & Mazin, D. 2009, A&A, 498, 25 

Salvaterra, R., & Ferrara, A. 2003, MNRAS, 339, 973 

Santos, M. R., Bromm, V., & Kamionkowski, M. 2002, MNRAS, 336, 1082 

Sullivan, L, et al. 2007, ApJ, 657, 37 

Thompson, R. L, Eisenstein, D., Fan, X., & Rieke, M. 2007a, ApJ, 657, 669 

Thompson, R. L, Eisenstein, D., Fan, X., & Rieke, M. 2007b, ApJ, 666, 658 

Tody, D. 1986, Proc. SPIE, 627, 733 

Tsumura, K., & Wada, T. 2011, PASJ, 63, 741 

Totani, T., & Yoshii, Y. 2000, ApJ, 540, 81 

Totani, T., et al. 2001, ApJ, 550, L137 

Vrtilek, J. M., & Hauser, M. G. 1995, ApJ, 455, 677 

Wada, T., et al. 2007, PASJ, 59, S515 



This manuscript was prepared with the A AS IATjtX macros v5.2. 



-24 - 



Table 1: Parameters used and obtained in this study 



Item 


2.4/iin 


3.2/zm 


4.1/iHl 


Number of stacked images 


40 


39 


24 


Total integration time (minutes) 


29.6 


28.9 


17.8 


Percentage of remaining pixels (%) 


47.7 


46.7 


47.1 


Limiting magnitude (AB mag) 


22.9 


23.2 


23.8 


Average sky brightness (nWm~ 2 sr _1 ) a 


114 


73 


105 


1 u fluctuation, sky (nWm _2 sr _1 ) a 


2.91 


1.72 


0.96 


1 a fluctuation, dark (nWm _2 sr~ 1 ) a 


1.98 


1.15 


0.59 



a Systematic errors are not taken into account. 



Table 2. Fluctuation of the observed sky, fluctuation due to faint galaxies (shot noise) and excess fluctuation in units 

of nWm _2 sr~ 1 



angle (arcsec) 




2.4 M m 






3.2 (im 






4.1 M m 




sky fluctuation 


shot noise 


excess fluctuation 


sky fluctuation 


shot noise 


excess fluctuation 


sky fluctuation 


shot noise 


excess fluctuation 


111 


0.35 ± 0.082 


0.13 ± 0.0026 


0.33 ± 0.089 


0.16 ± 0.035 


0.08 ± 0.0016 


0.14 ±.041 


0.078±0.02 


0.045±0.0009 


0.063±0.025 


125 


0.33 ± 0.043 


0.12 ± 0.0017 


0.31 ± 0.046 


0.16 ± 0.041 


0.075 ± 0.0012 


0.14 ± 0.047 


0.048 ± 0.011 


0.039 ± 0.0006 


0.027 ± 0.019 


143 


0.26 ± 0.035 


0.10 ± 0.0015 


0.24 ± 0.038 


0.13 ± 0.017 


0.064 ± 0.0010 


0.11 ± 0.020 


0.040 ± 0.0075 


0.034 ± 0.0005 


0.021 ± 0.014 


167 


0.20 ± 0.050 


0.087 ± 0.0021 


0.18 ± 0.056 


0.079 ± 0.034 


0.055 ± 0.0014 


0.057 ± 0.047 


0.039 ± 0.019 


0.030 ± 0.0007 


0.025 ± 0.029 


200 


0.093 ± 0.023 


0.070 ± 0.0012 


0.061± 0.035 


0.082 ± 0.016 


0.044 ± 0.0007 


0.068 ± 0.019 


0.031 ± 0.0070 


0.025 ± 0.0005 


0.019 ± 0.011 


250 


0.29 ± 0.045 


0.058 ± 0.0015 


0.29 ± 0.046 


0.13 ± 0.029 


0.036 ± 0.0009 


0.12 ± 0.030 


0.047 ± 0.010 


0.020 ± 0.0005 


0.042 ± 0.011 


333 


0.29 ± 0.12 


0.041 ± 0.0013 


0.29 ± 0.12 


0.13 ± 0.046 


0.027 ± 0.0009 


0.13 ± 0.047 


0.029 ± 0.012 


0.016 ± 0.0005 


0.025 ± 0.015 



Note. — 



Errors indicate 1 a error, i.e. standard deviation divided by square root of the number of data points for each ensemble. 
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Fig. 1. — Fluctuation maps obtained by stacking. The diameter of the circles is 600". Maps 
correspond to 2.4, 3.2, and 4.1 /^m from the left to the right. Upper panels indicate sky 
maps, while lower panels show dark maps. The unit is nWrn^sr -1 . The color scales shown 
in the bar below each map are chosen so that sky maps and dark maps cover the same ranges 
of sky brightness. 
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Fig. 2. — Amplitude maps of the two-dimensional fluctuation spectra, [g 2 P 2 (q)/(2vr)] 1 / 2 , in 
Fourier space. The results for the 2.4, 3.2, and 4.1 /im band are shown from the left to the 
right. The gray scales bars below each map indicate amplitudes of the fluctuation in units 
of nWm -2 sr _1 . 
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Fig. 3. — The upper panel shows the one dimensional fluctuation spectra, [q 2 P 2 /(2ii)} 1 / 2 in 
unit of nWm _2 sr _1 , obtained by two-dimensional Fourier analysis as a function of angular 
scale (2n)/q. Graphs correspond to the 2.4, 3.2, and 4.1 /j,m band from the left to the 
right. Filled circles and open triangles show the fluctuation spectra for sky and dark maps, 
respectively. The lower panel shows the fluctuation spectra of the sky after subtracting those 
of dark maps in quadrature. The straight lines indicate the fluctuation spectra of shot noise 
due to unresolved faint galaxies. All error bars represent 1 a error. 
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Fig. 4. — Results of subset analysis. The squares indicate the fluctuation spectra for the 
subtracted image for the case A subset, while the asterisks for the case B subset. The 
triangles show the fluctuation spectra for the dark maps. 
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Fig. 5. — Original fluctuation spectra (filled circles) are compared with those using common 
mask (open squares). The results for the 2.4, 3.2, and 4.1 fim band are shown from the left 
to the right. 
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Fig. 6. — The upper panel shows smoothed sky maps obtained by averaging pixels within 
a 50" diameter circle centered on each pixel. The lower panel shows smoothed dark maps 
obtained by the same procedure for the dark maps. The color scales shown in the bar 
below each map are chosen such that sky maps and dark maps have the same range of sky 
brightness. 
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Fig. 7. — The spectrum of the average fluctuation at large angles (100" ~ 350") (shown by 
filled circles) is compared with the Spitzer results (open squares). The open circles represent 
the spectrum of the correlating component normalized to the 2.4 /zm band. The solid line 
indicates a Ra y leigh- Jeans spectrum (~ A -3 ), while the dotted line indicates the spectrum 



in Figure 20 of 



Fernandez et al. 



(120101 ) . The vertical bars show 1 o error. 
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Fig. 8. — Two-point correlation functions, C 2 a{0), £3.2(0), and C,4.i(0), are shown from 
the top to the bottom. Filled circles and crosses correspond to the two-point correlation 
functions for sky and dark maps, respectively. Error bars indicate statistical errors, i.e. 
standard deviation divided by square root of number of pairs. 
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Fig. 9. — Point-to point correlation between wavelength bands for smoothed maps. The 
upper panel shows the correlation diagram for the 3.2 and 2.4 /xm bands, while the lower one 
shows that for the 4.1 and 2.4 /xm bands. Only independent data points are plotted. Error 
bars represent the standard deviations of the smoothed dark maps. The straight lines show 
the results of linear fits. 
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Fig. 10. — Cross correlation functions, C2.4(g>3.2(#) (upper panel) and (72.4<g>4.i(0) (lower panel). 
Error bars indicate statistical errors, i.e. standard deviation divided by square root of number 
of pairs. 
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Fig. 11. — The volume emissivity of zodiacal light at 2.2 (solid line) and zodiacal emission 
at 12 fxm (dashed line) along the line of sight towards the NEP. The vertical axis is normalized 
such that the total volume emissivity at the origin is 1.0. The top curve shows the total 
emissivity that is, sum of the smooth cloud, the dust bands, and the circumsolar ring. The 
middle and bottom curves indicate the volume emissivities of the circumsolar ring and the 
dust bands. 
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Fig. 12. — Diagrams of correlations between AKARI maps and the far-infrared map (90 fim) 
observed with FIS/AKARI. The angular resolutions of the AKARI maps were degraded to 
that of the FIR observations, s~ 30". Error bars indicate the standard deviations of the 
dark maps for both observations. 



